%%% PREPROCESSING  %%%

close all
clear all

% directories - adjust depending on your own data organization
rootdir = '/home/predatt/anatod/ISI'; % root directory
datadir = '/home/predatt/anatod/ISI/raw';  % raw data directory
gavgdir = '/home/predatt/anatod/ISI/MEG/GAVG'; % grand average directory

% load subject and dataset information
cd(gavgdir)
load subjname
load dataset

clc

% enter subjects to be processed
subjects = 1;

for sb = subjects;

%% define the trials
infofile = strcat(subjname{sb},'_datainfo.mat'); % create file to fill with participant information
subdir = fullfile(rootdir,'MEG',subjname{sb}); % subject directory
fprintf('Processing subject %s... \n', subjname{sb});

cd(datadir) % 

    cfg = [];
    cfg.dataset = dataset{sb};
    cfg.trialdef.eventtype  = 'UPPT001';
    cfg.trialdef.eventvalue = [1:8 17:24 33:40 49:56 65:72]; % read all conditions at once
    cfg.trialdef.prestim    = 0.5;
    cfg.trialdef.poststim   = 1;
    cfg.trialfun = 'ft_trialfun_general';
    
    cfg.dataformat= 'ctf_meg4';
    cfg.headerformat= 'ctf_res4';
    
    [cfg] = ft_definetrial(cfg);
    
%% additional trial information
    
    % add tone2 triggers from each trial
    cfg.trl(:,5) = cfg.trl(:,4)+8;

    % add reaction times
   event = cfg.event;
   sel = ismember({event.type}, 'UPPT001');
   valchar    = ischar(cfg.trialdef.eventvalue); 
   valnumeric = isnumeric(cfg.trialdef.eventvalue); 
  for i=1:numel(event)
    if (ischar(event(i).value) && valchar) || (isnumeric(event(i).value) && valnumeric)
      sel(i) = sel(i) & ~isempty(intersect(event(i).value, cfg.trialdef.eventvalue));
    end
  end
  
sel = find(sel);
smprate = 1200;
resp_trig = 1; % response button
% resp_trig = 201; % response button is 201 for subjects 21-24

rtdata = [];

    % find all responses within 3 seconds from stim1
for i=sel
        next_val = [event(i+1:end).value];
        for j=1:length(next_val)
            if next_val(j) == resp_trig; 
                if event(i+j).sample - event(i).sample < 3*smprate;
                    resp = 1;
                    RT = (event(i+j).sample - event(i).sample)/smprate;
                    break;
                else resp = 0; RT = NaN;
                    break;
                end
            else resp = 0; RT = NaN;
            end
        end
    
    rtdata = [rtdata; [resp RT]];

end



cfg.trl(:,6:7) = rtdata;
    
    
%% create the datainfo structure with all relevant information
    
        datainfo = {};
        datainfo.basicinfo.dataset = cfg.dataset;
        datainfo.basicinfo.datafile = cfg.datafile;
        datainfo.basicinfo.headerfile = cfg.headerfile;
        datainfo.basicinfo.dataformat = cfg.dataformat;
        datainfo.basicinfo.headerformat = cfg.headerformat;
        datainfo.trialfun = cfg.trialfun;
        datainfo.trialdef = cfg.trialdef;
      
        
    %% calculate reaction times and accuracy
    
    deviant_trigs = intersect(cfg.trialdef.eventvalue, [cfg.trialdef.eventvalue] + 4);
    rt_deviants = cfg.trl(find(ismember(cfg.trl(:,4),deviant_trigs)),:);
    hits = sum(rt_deviants(:,6)); % how many hits?
        
    standard_trigs = intersect(cfg.trialdef.eventvalue, [cfg.trialdef.eventvalue] -4);
    rt_standards = cfg.trl(find(ismember(cfg.trl(:,4),standard_trigs)),:);
    falsealarms = sum(rt_standards(:,6)); % how many false alarms?
        
    att1_weakjitter_trigs = [5 21 37 53 69]; % condition 1
    att2_weakjitter_trigs = [5 21 37 53 69]+1; % condition 2
    att1_strongjitter_trigs = [5 21 37 53 69]+2; % condition 3
    att2_strongjitter_trigs = [5 21 37 53 69]+3; % condition 4

    rt_deviants_cond1 = rt_deviants(find(ismember(rt_deviants(:,4),att1_weakjitter_trigs)),:);
    mean_rt_cond1 = nanmean(rt_deviants_cond1(:,7));

    rt_deviants_cond2 = rt_deviants(find(ismember(rt_deviants(:,4),att2_weakjitter_trigs)),:);
    mean_rt_cond2 = nanmean(rt_deviants_cond2(:,7)) - 0.5;

    rt_deviants_cond3 = rt_deviants(find(ismember(rt_deviants(:,4),att1_strongjitter_trigs)),:);
    mean_rt_cond3 = nanmean(rt_deviants_cond3(:,7));

    rt_deviants_cond4 = rt_deviants(find(ismember(rt_deviants(:,4),att2_strongjitter_trigs)),:);
    mean_rt_cond4 = nanmean(rt_deviants_cond4(:,7)) - 0.5;
    
    % add rt information to datainfo
        datainfo.rt.conditions = {'att1 exp 500 ms', 'att 2 exp 500 ms', 'att1 exp anytime', 'att2 exp anytime'};
        datainfo.rt.rt = [mean_rt_cond1, mean_rt_cond2, mean_rt_cond3, mean_rt_cond4];
        datainfo.rt.hits = hits;
        datainfo.rt.falsealarms = falsealarms;
        datainfo.rt.numstandards = 486; % total number of standards
        datainfo.rt.numdeviants = 48; % total number of deviants
    
    %% remove deviant trials and trials with responses
    
    cfg.trl = rt_standards;
    ind_falsealarms = find(rt_standards(:,6)==1);
    cfg.trl(ind_falsealarms, :) = [];
    cfg.trl = cfg.trl(:,1:5);
    cfg.trlorig = cfg.trl;
    
    datainfo.trlorig = cfg.trl;
    datainfo.trl = cfg.trl;
   
    %% save, save, save!

    
 cd(subdir)
 save(infofile,'datainfo');
 close all
 
 % the datainfo structure is later used to define which parts of the raw data to load

end


